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Abstract 

Simultaneous localization and tracking (SLAT) in sensor networks aims to determine the positions 
of sensor nodes and a moving target in a network, given incomplete and inaccurate range measurements 
between the target and each of the sensors. One of the established methods for achieving this is 
to iteratively maximize a likelihood function (ML), which requires initialization with an approximate 
solution to avoid convergence towards local extrema. This paper develops methods for handling both 
Gaussian and Laplacian noise, the latter modeling the presence of outliers in some practical ranging 
systems that adversely affect the performance of localization algorithms designed for Gaussian noise. 
A modified Euclidean Distance Matrix (EDM) completion problem is solved for a block of target 
range measurements to approximately set up initial sensor/target positions, and the likelihood function is 
then iteratively refined through Majorization-Minimization (MM). To avoid the computational burden of 
repeatedly solving increasingly large EDM problems in time-recursive operation an incremental scheme is 
exploited whereby a new target/node position is estimated from previously available node/target locations 
to set up the iterative ML initial point for the full spatial configuration. The above methods are first 
derived under Gaussian noise assumptions, and modifications for Laplacian noise are then considered. 
Analytically, the main challenges to be overcome in the Laplacian case stem from the non-differentiability 
of l\ norms that arise in the various cost functions. Simulation results confirm that the proposed algorithms 
significantly outperform existing methods for SLAT in the presence of outliers, while offering comparable 
performance for Gaussian noise. 
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I. Introduction 

This work addresses the problem of tracking a single target from distance-like measurements taken by 
nodes in a sensor network whose positions are not precisely known. The goal is to estimate the positions 
of all sensors and the target, given only partial or no a priori information on the spatial configuration 
of the network. The ability to track a target is a key component in several scenarios of wireless sensor 
networks, hence methods that avoid the need for careful calibration of sensor positions are practically 
relevant. 

In 12, (21 SLAT is formulated in a Bayesian framework that resembles the related and well-studied 
problem of Simultaneous Localization and Mapping (SLAM) in robotics. The a posteriori probability 
density function of sensor/target positions and calibration parameters is recursively propagated in time as 
more target sightings become available. In [1] these observations are true range measurements obtained 
through a combination of transmitted acoustic and radio pulses. Alternatives to range include pseudorange 
and bearing information estimated from camera images O, or the (somewhat unreliable) Received 
Signal Strength (RSS) of radio transmissions 0. In (4[ the SLAT problem is also formulated in a 
Bayesian framework as a general state evolution model under a binary proximity model and solved 
in a decentralized way using binary sensor networks. Other SLAT-like approaches include localization 
(calibration) as presented in Q, where positions and orientations of unknown sources and sensors are 
centrally obtained via ML based on time-of-arrival and angle-of-arrival measurements. Several of the 
above references emphasize decentralized processing, in line with the local observation model. 

When target dynamics are not accounted for, the SLAT problem may be thought of as a special case 
of Sensor Network Localization (SNL) with a limited set of intersensor measurements. In fact, EDM 
and related matrix completion methods based on Semidefinite Programming (SDP) have been adopted 
previously for static SNL (see [6] and references therein). EDM completion for SLAT has been discussed 
in (31, although the authors ended up pursuing an alternative approximate completion approach based on 
a variant of Multidimensional Scaling (MDS). Underwater and underground scenarios with uncertainty 
in anchor positions are considered in Q, and edge-based SDP is proposed to reduce the computational 
complexity of SNL. In [ 8 ] static SNL is formulated as a problem of ML phase retrieval. 

In addition to centralized SNL approaches such as O-Hl, enumerated above, a wealth of results is 
available on distributed approaches for scenarios where the existence of a central node is inconvenient, 
e.g., due to congested communications in its vicinity, or excessive vulnerability of the whole infrastructure 
to failure of that single node (3J, |[9l- lfl3l . A two-step approach based on second-order cone programming 
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relaxation with inaccurate anchor positions is introduced in Q- In ifTTl a weighted least-squares algorithm 
with successive refinement provides both position estimates and their covariances in partially connected 
scenarios. A distributed weighted MDS method with majorization approximations is applied in 1121 . The 
cost function and majorization technique are similar to the ones used in this paper for ML iterative 
refinement under Gaussian noise, but initialization relies on prior estimates of sensor positions. A 
wholly different iterative approach for distributed SNL using barycentric coordinates and Cayley-Menger 
determinants is developed in |[T3l . For m -dimensional Euclidean space a node reinterpolates its position 
based on estimates from a set of m + 1 neighbors such that it lies in their convex hull. 

This paper focuses on centralized SLAT based on plain ML estimation. A basic iterative optimization 
approach using Majorization Minimization (MM) is first developed for batch estimation, i.e., when all 
measurements are processed simultaneously. A time recursive method is then obtained by estimating 
each new target or sensor position as the corresponding range measurements become available, and 
then iteratively re-optimizing the expanded ML cost function. This continuation scheme only requires 
initialization of all target/sensor positions at the first time step, which is computationally less complex 
than doing so for each new sensed target or sensor position. EDM completion for batch estimation is 
proposed to initialize the iterative ML algorithm with little a priori knowledge of sensor/target positions, 
thus alleviating the problem of convergence to undesirable local extrema. 

This approach was proposed in lfl4l for Gaussian noise, using cost functions for batch and time 
recursive initialization schemes that match squared observations with squared estimated ranges. These 
discrepancies are eliminated in the present paper, such that both initialization and ML refinement operate 
with cost functions that match plain (non-squared) ranges, leading to improved robustness under strong 
measurement noise lfl31l . We use the Source Localization in Complex Plane (SLCP) approach proposed in 
Ifl6l to obtain a cost function for initialization of incremental target/sensor position estimates that admits 
an accurate convex relaxation as an SDP. 

This work also develops modified versions of initialization and ML refinement for Laplacian noise, 
which models the presence of outliers in some practical ranging systems that adversely affect the 
performance of localization algorithms designed for Gaussian noise HI, IfTTl . This is accomplished by 
replacing £2 norms with l\ norms in the optimization problems for initialization and refinement that were 
previously formulated for Gaussian noise, and then performing suitable manipulations to write these in a 
form that is amenable to general-purpose solvers. A source localization algorithm with i\ norms (SL^i) 
is also derived for incremental initialization of the SLAT scheme. 

In IfTTl i\ norms are also used to handle outliers, but the proposed method is very different from the 
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one developed here, as it relies on linear programming to identify the outliers, and then remove them 
from consideration when computing the source location. In our work all measurements are kept, as the 
modulus of range differences that appears in cost functions ensures that outlier terms do not overwhelm 
the remaining ones as long as the proportion of outliers remains small. Another approach for handling 
outliers is presented in Ifl8l , where the Huber cost function interpolates between i\ and £2 norms. This 
function is minimized via iterative majorization techniques with a priori information on sensor positions, 
where in each step a majorization subproblem is solved using Costa's algorithm. 

The paper is organized as follows. In section [TTJ the SLAT problem is introduced. Section UlTI presents 
estimation methods for range measurements corrupted by Gaussian noise, namely, EDM initialization, 
iterative likelihood refinement by MM, and time -recursive updating through incremental estimation of 
target/sensor positions. Section [TV] develops similar methods for Laplacian noise. Section [V] provides 
simulation results of batch and time recursive approaches. In addition, the performances of initialization 
techniques are compared with and without outliers. Finally, section [VTI summarizes the main conclusions. 

Throughout, both scalars and individual position vectors in (2D) space will be represented by lowercase 
letters. Matrices and vectors of concatenated coordinates will be denoted by boldface uppercase and 
lowercase letters, respectively. The superscript T (H) stands for the transpose (Hermitian) of the given 
real (complex) vector or matrix. Below, I rn is the mxm identity matrix and l m is the vector of m ones. 
For symmetric matrix X, X y means that X is positive semidefmite. 



The network comprises sensors at unknown positions {xx, x%, • • • , x n } G IR 2 , a set of reference sensors 
(anchors) at known positions {ai, 02, • • • , a«} G K 2 , and unknown target positions {ei, e2, . . . , e m } G R 2 . 
A central processing node has access to range measurements between target positions and all sen- 
sors/anchors, namely, 



where Wij and Wkj denote noise terms. A practical system that provides such range measurements is 
used, e.g., in UJ. 

a) SLAT Under Gaussian Noise: If disturbances are Gaussian, independent and identically dis- 
tributed (i.i.d.), maximizing the likelihood for the full batch of observations is equivalent to minimizing 
the cost function 



II. Problem Formulation 
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The full set of unknown sensor and target positions is concatenated into column vector x € R 2 ( n+m ), the 
argument of Q,g- The goal of our SLAT approach is to find the set of coordinates in x which minimizes 



b) SLAT Under Laplacian Noise: When the disturbances are Laplacian and i.i.d., thus heavier tailed 
than Gaussian, maximizing the likelihood amounts to minimizing the cost function 



When compared with (HJ, the absence of squares in the summation terms of © renders the function less 
sensitive to outlier measurements dij with large deviations from the true ranges. 

Due to the nature of this problem the functions Qg an d are invariant to global rotation, translation 
and reflection in the absence of anchors. In order to remove ambiguities in the solutions, at least I = 3 
non collinear anchors must be considered. As in many other ML problems, the functions &g and are 
in general nonconvex and multimodal, hence their (approximate) minimization proceeds in two phases, 
initialization and refinement. The former provides a suitable initial point through EDM completion, which 
tends to avoid convergence towards undesirable local minimizers of the ensuing iterative refinement 
algorithms based on MM or weighted-MM. 



This section develops algorithms for EDM initialization, MM refinement, and recursive estimation 
in SLAT under the assumption that measurement noise is i.i.d. and Gaussian. Modifications for i.i.d. 
Laplacian noise are considered in Section [TV] A basic formulation of EDM completion with squared 
distances is provided first, and will form the basis for the initialization methods described in Sections 
HIFBl and HV-Al 

A. EDM with Squared Distances 

The basic EDM completion problem, described below, operates on squared ranges lfl9l , ll20l . Even 
though it is not matched to the likelihood function £[]), it will be useful for benchmarking in Section |VJ 
as its performance is representative of other popular SNL methods 0, II2T1 and the SLAT approach of 



A partial pre-distance matrix D is a matrix with zero diagonal entries and with certain elements 
fixed to given nonnegative values which are the squared observed distances, = d?j. The remaining 
elements are considered free. The nearest EDM problem is to find an EDM E that is nearest in the 



©• 




(2) 



III. SLAT under Gaussian Noise 
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least-squares sense to matrix D, when the free variables are not considered and the elements of E satisfy 
Eij = \\yi — yj\\ 2 for a set of points yi. The geometry and properties of EDM (a convex cone) have been 
extensively studied in the literature Ifl9l , EOl . The nearest EDM problem in 2D space is formulated as 

minimize £ ije0 {Eij - d\ ) 2 
E 

subject to E G £ (3) 
E(A) = A 
rank(JEJ) = 2, 

where 

J = (l p - il p l p T ), p = m + n + l, 

is a centering operator which subtracts the mean of a vector from each of its components. In ©, O is the 
index set for which range measurements are available. The constraint E(^4) = A, where A is the index 
set of anchor/anchor distances and A^ = ||aj — dj\\ 2 is the corresponding EDM submatrix, enforces the 
known a priori spatial information. Matrix E belongs to the EDM cone £ if it satisfies the properties 

En = 0, Eij > 0, -JEJ y 0. (4) 

The rank constraint in © ensures that the solution is compatible with a constellation of sensor/anchor/target 
points in R 2 . Extraction of the set yi from E is described below. Problem © is also known as the penalty 
function approximation |[T9l due to the form of the cost function v?i(E) = Ylij(Eij ~ dfj) 2 - Expressing 
© in terms of full matrices and dropping the rank constraint, a compact relaxed SDP formulation is 
obtained as 

minimize ||W©(E-D)||| 

E (5) 

subject to E G £, B(A) = A, 
where W is a mask matrix with zeros in the entries corresponding to free elements of = d? , and 
ones elsewhere. Combined with the Hadamard product 0, and Frobenius norm this replaces the 

summation in © over the observed index set O. From here on, we will call this method EDM with 
squared ranges (EDM-SR). 
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B. SLAT Initialization: EDM with Plain Distances 

Instead of trying to match squared distances, we can apply EDM completion to plain distances as 

minimize YlijiV^ij ~ d ij) 2 



E 

(6) 

subject to E G £, B(A) = A 
rank(JEJ) = 2. 

For this method the penalty function is ^2(13) = £\ . ( E^j — d^) 2 , which more closely matches the 



terms in the likelihood function (fl}, and © is thus expected to inherit some of the robustness properties 
of ML estimation. Expanding the objective function in Q results in 

minimize J2i,j( E ij ~ ^\/E~dij + d? ) 
E 

(V) 

subject to E e £, B(A) = A 
rank(JEJ) = 2. 

A relaxed SDP is obtained by introducing an epigraph-like variable T and dropping the rank constraint 

minimize J2i,j(. E ij ~ 2T ijdij) 
E,T 



subject to Tf- < Ei 



(8) 

E G £ , E(A) = A. 
From here on, we will call this method EDM with plain ranges (EDM-R). 

Remark that the solutions of the initialization techniques described here and in Sections [III-AI and IIV-AI 
are distance matrices. Detailed explanations of how to estimate the spatial coordinates of the sensors and 
target positions from EDM and the usage of anchors are given in lfl4l . The basic idea is to use a linear 
transformation to obtain the Gram matrix Y T Y from the EDM matrix E, from which spatial coordinates 
Y are extracted by singular value decomposition (SVD) up to a unitary matrix. The anchors are then used 
to estimate the residual unitary matrix by solving a Procrustes problem. As discussed in lfl4ll . observation 
noise can significantly disrupt the estimated sensor/target coordinates through EDM completion and rank 
truncation, and it was found that much more accurate results are obtained by using those as a starting point 
for likelihood maximization. MM algorithms are proposed next for iterative likelihood maximization. 

C. SLAT Iterative Refinement: Majorization-Minimization 

The key idea of MM is to find, at a certain point x*, a simpler function that has the same function 
value at x* and anywhere else is larger than or equal to the objective function to be minimized. Such a 
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function is called a majorization function. By minimizing the majorization function we obtain the next 
point of the algorithm, which also decreases the cost function [22 J. The detailed derivation of MM is 
given in [14], but it is summarized below for the reader's convenience. Define two convex functions as 

fij(x) = \\xi - ej\\, 9kj(x) = \\a k - ej\\. (9) 

Expanding / and g in £T|) and using first-order conditions on convexity |fT9l . 

n G (x) < - to xi ) + < v ^( x *)> ( x - x *)» + 4) 



(10) 



+ J2(9kj(x) ~ 2c4i(2fcj(x i ) + (V^x 4 ), (x - x*))) + d%) , 

where (u, v) = u T v, we get the proposed majorization function on the right side of ([TOl . which is 
quadratic in x and easily minimized. Hence the MM iteration 



x t+1 = argmin^r (/g(x) - 2d ii (V/y(x*) J x)) + £ (<^(x) - 2d kj (Vg kj (x t ), x)) (11) 
turns out to be obtained as the solution of a linear system of equations. 



D. Time-Recursive Position Estimation: SLCP 

Suppose that a batch of observations has been processed, and a new target position is to be estimated. 
We could repeat the previous approach by redefining the new batch as the old one concatenated with 
the new set of observations. However, this is computationally expensive due to the EDM step. Also, 
previous estimated positions would be ignored. To alleviate the load we propose a simple methodology 
to obtain a good initial point which avoids the EDM step. It consists of fixing the previous positions at 
their estimated values and only estimating the new target position. More precisely, we minimize 

n+l 

=1^(11^-2/11 -^f, (12) 
i=l 

where y is the new target position, b\ denotes the position of a sensor or anchor, and d{ is the corresponding 
range measurement. We propose SLCP [16] to minimize (TT2l . and briefly summarize the method below. 

One can readily see that each term in ([L2l quantifies the distance between the target location and a circle, 
centered at an anchor or sensor, with radius equal to the measured range, e.g., | \\bi — y\\ — dA = \\y — yi\\, 
where yi is located at the intersection of the line connecting fej and y with the circle {z : \\z — bi\\ = di}. 
An equivalent formulation is therefore 

En+l II n 2 

i=i \\y - ViW 

y, yi (13) 

subject to \\bi — y%\\ = di i = 1, . . . , n + I. 
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SLCP relies on a complex formulation of (fT3l ) for the 2D case, and the problem is then manipulated and 
relaxed to an SDP of the form 

maximize t H — Vrr T <I>r 

(14) 

subject to * >z 0, 4>u = 1 
Ac H $c > t 2 , 

where r and c are constant vectors built from sensor/anchor complex coordinates and measured radii. 
The solution of the optimization problem ([141) is a positive semidefinite matrix, hopefully with near-1 
rank. Afterwards, the target coordinates are estimated by SVD of $ as described in lfl4l . 

After an optimal target position is obtained, we return to the cost function £[]) or © and iteratively 
refine all the estimates. This incremental or time recursive procedure can be applied to either new targets 
or sensors. 

In a previous paper lfT4l . the source localization method derived in |[23l , termed Squared Range Least 
Squares (SR-LS), is proposed as the time recursive method. Note that in |[23l squared distances are 
matched, leading to a Trust Region optimization problem. However, as demonstrated in fTo*! . SLCP is a 
more accurate source localization method and its cost function (fT2l) is better matched to the likelihood 
function £T|). This makes it more convenient for initialization of iterative refinement algorithms, which 
therefore require fewer iterations to converge and/or are less likely to be trapped in undesirable local 
extrema. 



IV. SLAT under Laplacian Noise 

A. SLAT Initialization: EDM with Ranges and l\-norm 

Among the penalty function approximation methods ^i-norm approximation is known to be robust 
to outliers fl9l . Therefore, the penalty function of the third SLAT initialization method is chosen as 
(/33(E) = j I \J Eij — d%j\, and the associated optimization problem becomes 

minimize £\ ^ | ^/E^ - dij \ 
E 

(15) 

subject to E G £ , E(A) = A 
rank(JEJ) = 2. 

This problem is not convex, as the objective function | ^jEij — dij \ is convex when \JTkj — dij < 0, but 
concave for J Eij — d^ > 0. To obtain a convex approximation the objective function is replaced by a 
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Fig. 1 : The value of | y/Eij — d™ | vs E^ , and the linear approximation of the concave part. 



linear approximation 

&ijEij -\- bij , 



4 



\Z-^max ~\~ d{j y/E max ~i~ dij 

in part of the domain where it is concave, as shown in Figure [TJ The two functions coincide for E± 



(16) 



d 2 - 



(17) 



and Eij = E max , where the constant E max is a practical upper bound on (squared) range measurements. 

I E~ij , QijEij+bij} and use the epigraph 



Thus, we replace | y/ E^ — d^ | by its convex envelope maxjdjj - 
variable T to obtain 

minimize YU,j T v 
E, T 

subject to maxjdjj — ^Eij,aijEij + bij} < Ty 
E € £ , B(A) = A 
rank(JEJ) = 2. 
A relaxation of (fTTT ) after dropping the rank constraint is 

minimize Yi,j T ij 
E, T 

subject to {d^ — Tij) 2 < Ejj 
Q'ijEij ~t~ bij ^ Tij 
E e £, B(A) = A. 



(18) 



Note that the first constraint in (1181 ) is not equivalent to d 



'Eij < Tij, but rather to 



1 En < 



dij — T^ < yj Eij , which amounts to intersecting the original epigraph with the parabolic hypograph 
dij + \jEij > T^ . This preserves the convexity of the feasible set and does not change its lower boundary 
for E^ G [0, E max ], where the optimal point will be found. The constraint can now be readily expressed in 
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standard form without introducing additional variables, e.g., as an LMI or a second-order cone constraint 



i d. L j T^j 


t o 




2(d ij 


-Tij) 


or 








Eij 


- 1 



This technique will be called EDM with ranges and ^-norm (EDM-R-£i). 



< Eij + 1. 



(19) 



B. SLAT Iterative Refinement: Weighted Majorization Minimization 

Robustness to outliers in the cost function (fSJ) for Laplacian noise is gained at the expense of differen- 
tiability. To circumvent that shortcoming we resort to the well-known re-weighted least squares approach 
ll25l . which replaces the minimization of (f2]) with a sequence of minimizations of smooth approximation 
functions that converge to Ql(x). Specifically, © is first written as 

II j \2 

1 1 1 ) i 



n L (x) = Y, 



Uij 



ej\\ ~ dij) 2 + ^v k j( 



\a k 



(20) 



with 



1 



di 



1 



dkj 



At time t the function to be minimized becomes f2^(x), which has the same form of d20l but the functions 
Uij , Vkj above are now replaced by constants based on the estimated positions after the previous iteration 

1 , 1 



d; 



dkj 



(21) 



An inner optimization loop could now be used to minimize (x) for every t but, as shown in Appendix 
lAl a single iteration suffices to ensure convergence. With fixed u\j, u| . the same majorization technique 
of Section IIII-CI then yields the weighted-MM iteration 

x m = argmin^4 (/?.(x) - 2d ij (V^(x*), x)) +J2 v kj (^'( x ) " 2d fci (V 5 fci(x t ) ) x)) . (22) 
hi k,j 

In practice the weights u\j and v\j must be modified to avoid the possibility of division by zero |[26l . 
which in our case is achieved by saturating them at 10 5 when computing (l22l . 



C. Time-Recursive Position Estimation: SL£\ 

The ML source localization problem under Laplacian noise is equivalent to 

minimize ^ L {y) = Ya=i 1 1 1 V ~ h 1 1 ~ <U \ 

y 



(23) 
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or 

minimize (£^=1 I \\v - bill - di\) 2 , 

y 

where y, 6, and d{ are denned in section ITlI-DI We use ideas from |27] to express the minimization of 
^ 2 L as a weighted sum of squares. 



Lemma IV-C.l. The following problem is equivalent to (124b 

■ ■ ■ ■ ■ ■ ^n+l (\\y-bA\-diY 

minimize minimize 2^j=i — \. 

y A g M n+i (25) 
subject to Aj > 0, 1 T A = 1. 

A proof is given in Appendix |B] As claimed in section IIII-DI the difference between the true range 
and observed range is actually equivalent to the distance between the source position and the point on 
the circle with center bi and radius di. An equivalent formulation is therefore 

En+l \\y— Vi\\ 2 
i=l x 

V ^ X (26) 
subject to \\yi — bi\\ = d{ 

Xi > 0, 1 T X = 1. 

If we fix the yi and Xi, the solution of (1261) with respect to y is an unconstrained optimization problem 
whose solution is readily obtained by invoking the optimality conditions 

n+l , , v-^n+Z yi_ 

Geometrically, the first constraint of (1261) defines circle equations, which can be compactly described 
in the complex plane as y; t = + d{e^\ We collect these into a vector y = b + Ru, where b = 



h ... b n+ i 



T 

£ C n+l , R = diag(di,..., d n+ i) £ r(«+0x(«+0, and u 
C n+i . Using the optimal y, we get 

minimize y^IIy = (b + Ru) H II(b + Ru) 
A, u 



r 



( 28 ) 

subject to A, > 0, 1 T A = 1 

I it,: I = 1, 
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where 



n 








l 

A„ +i 



1 

En+l 1 



_1_ 



1 

A„ +i 



L Ai ••• A n+! 



(29) 





'•■ 



= A- 1 -A- 1 l(l r A- 1 l)- 1 l T A- 1 > 
with A = diag(Ai, . . . , X n +i)- 

Matrix II resembles an orthogonal projector. Using the matrix inversion lemmaj it is seen to be the 
limiting case II = lim -_ i . 00 (A + all T )~ l and thus positive semidefinite. This format is more amenable 
to analytic manipulations in optimization problems and will be used throughout. The parameter a is taken 
as a sufficiently large constant (see Appendix |B]), although it could also be regarded as an additional 
optimization variable to ensure adequate approximation accuracy. 

We now introduce an epigraph variable t G K in (|28T ). i.e., we minimize over t and add the constraint 
t — (b + Ru)^n(b + Ru) > 0. Applying Schur complements the constraint may be successively written 
as 



>- 



! (b + Ru)(b + Ru 

t 



y o. 



(30) 



t ib • Hllj" -K TJ„WK ; t>„,// 

b + Ru n 1 

The formulation becomes 

minimize t 
A, u, t 

subject to Aj > 0, 1 T A = 1 

\Ui\ = 1 

tA + ta\\ T y (b + Ru)(b + Hu) H . 

Finally, we define B = [b R], v H = [1 u H ], V = vv , and ignore the rank- 1 constraint on the new 
variable V to obtain the relaxed SDP 

minimize t 



(31) 



p,v,t 

subject to fa > 0, 1 T (3 = t 

diag(/3) + tall T y BVB F . 



(32) 



1 (A + BCD)" 1 = A' 1 - A- 1 B{DA- 1 B + C~ 1 )- 1 DA- 1 . 
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The solution of the optimization problem (l32l includes the positive semidefinite matrix V from whose 



first row or column u can be extracted directljo to obtain y = b + Ru and the target coordinates from 



As in section IIII-DI after an optimal target position is obtained we return to the cost function (O and 
iteratively refine all the estimates. It is demonstrated in simulation in section [V] that SL^i is more robust 
to outliers than the SLCP algorithm of section Mi-Pi as its cost function (l23l is better matched to the 
likelihood function Q. 



Example 1 [Comparison of Batch Initialization Methods]: To investigate the accuracy of the 
methods, we set a physical scenario containing four anchors, five unknown sensors, and six target positions 
in a [0, 2] x [0, 2] area. Range measurements are corrupted by additive spatio-temporally white noise with 
standard deviation 

Ogaussian £ [0.005, 0.03]. This noisy observation model may lead to near-zero or negative 
range measurements, in which case we follow normal practice [6] and set them equal to a small positive 
constant (10~ 5 in our simulations). With the chosen noise variances this occurs sufficiently seldom (up 
to 0.04% of measurements) for its impact on estimation accuracy to be unimportant. Several algorithms 
are tested (EDM-SR, EDM-R, EDM-R-Zi, MM initialized by EDM-SR (EDM-SR+MM), MM initialized 
by EDM-R (EDM-R+MM) and MM initialized by EDM-R-Zi (EDM-R-Zi+MM)), and their performances 
compared according to the total root mean square error (RMSE) 



where x\ denotes the i-th estimated sensor or target position in the A;-th Monte Carlo run for the specific 
noise realization. In each of K = 150 Monte Carlo runs, a random network is generated according to 
the physical scenario described above. To assess the fundamental hardness of position estimation, error 
plots for Gaussian noise also show the total Cramer-Rao Lower Bound (CRLB), calculated as 



for each noise variance, where CRLB^ denotes the /c-th diagonal element of the matrix lower bound. 
The CRLB for anchored and anchor-free localization using ranging information has been studied in 

2 Alternatively, u can be obtained by rank 1 factorization of the lower right submatrix of V corresponding to uu ff , as in |6|, 




V. Numerical Results 




(33) 




(34) 



|21]. 
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I28l - ll30l for different variance models of range estimation noise. For convenience, the CRLB for our 
SLAT problem under Gaussian noise is rederived in Appendix O in terms of the notation adopted in 
this paper. We do not prove unbiasedness of our estimators, a mathematically challenging endeavor that 
would be required to fully justify benchmarking against the CRLB. In our experimental results, however, 
we found no clear evidence of bias for small noise levels, where convergence to undesirable extrema 
of the cost functions is avoided. Figure [2a] shows that plain EDM-R has better accuracy than EDM-SR 
and EDM-R-Zi, although the performance gap closes after iterative refinement by MM. Moreover, MM 
initialized by the various methods nearly touches the CRLB except when the noise variance is large. 

To compare the total RMSE of the algorithms in the presence of outliers, modified range measurements 
are created according to a "selective Gaussian" model di = || • || +Wi + |ej|, where e« is a white Gaussian 
noise term with standard deviation o" out ii e r £ [0.4,2]. The disturbance e« randomly affects only two range 
measurements, whereas Wi with <r gauss j an = 0.01 is present in all observations. This outlier generating 
model deviates from the earlier Laplacian asssumption, but it is arguably representative of observed range 
measurements in practical systems [8]. Numerical results under a pure Laplacian model will be presented 
in Examples 3 and 4. In the presence of high noise and/or outliers, Fig. [2b] shows that weighted-MM 
refinement does not close the performance gap between EDM-R-^i, EDM-R and EDM-SR initialization 
because in the latter cases the algorithms converge more often to local minima, thus producing a larger 
total RMSE. 

Example 2 [Uncertainty Ellipsoids]: To further examine the accuracy of MM and weighted-MM 
with different initialization methods, we randomly generated two networks of 10 sensors, 4 anchors and 
11 target positions. 100 Monte Carlo runs were used to find the mean and (la) uncertainity ellipsoids of 
the positions estimated by the methods. The mean and uncertainity ellipsoids for <7 gauss j an = 0.025 and 
fgaussian = 0.02/<7 out ij er = 0.8 are shown in Figs. [3] and [4] respectively. Again, outliers are randomly added 
to two range measurements in Fig. |4] 

Without outliers (Fig. using EDM-SR, EDM-R, or EDM-R-/i as an initialization to MM makes the 
uncertainity ellipsoids shrink dramatically after refinement, yielding very similar means and covariances. 
These are only displayed in the detail view of Fig. [3b] as they are too small to be shown in Fig. [3a] In 
the presence of outliers, (Fig. [4]) the uncertainity ellipsoids of EDM-SR+wMM are bigger than for other 
methods and the means of the estimated positions are shifted. Since EDM-R-Zi and EDM-R initializations 
converge to global extrema most of the time, the means of the positions estimated by weighted MM still 
approach the true positions and their uncertainity ellipsoids are much smaller than for EDM-SR+wMM. 
In the presence of outliers this example shows that EDM-R-Ii+wMM is clearly superior to the other 
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Fig. 2: Comparison of batch initialization and refinement methods for SLAT. 
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methods. 

A Note on Practical Computational Complexity of EDM Initialization: Our experiments were con- 
ducted on a machine with an Intel Xeon 2.93 GHz Quad-Core CPU and 8 GB of RAM, using Matlab 7.1, 
CVX 1.2 and Yalmip 3/SeDuMi 1.1 as a general-purpose SDP solver. CPU times are similar for EDM- 
SR, EDM-R and EDM-R-11, under 5 seconds for the example described above with n = 25 unknown 
positions and empirically increasing with n 4 5 for larger values of n (< 100). This gives a notion of what 
network sizes are currently practical for the EDM initialization methods, while keeping in mind that CPU 
times are known to be unreliable surrogates for intrinsic computational complexity due to dependencies 
on factors such as machine hardware architecture, operating system, efficiency of numerical libraries, 
and solver preprocessing. No attempt was made to formulate the EDM completion problems in the most 
efficient way possible for the SDP solver. For MM type iterative algorithms extremely large problem 
sizes can be efficiently handled using contemporary numerical algorithms and computing platforms. In 
our experiments each iteration takes up to about 1 millisecond. 

Example 3 [Comparison of Source Localization Methods for Time-Recursive Initialization ]: In 
this example SL^i and SLCP are compared using five anchors. We performed 100 Monte Carlo runs, 
where in each run the anchor and source positions were randomly generated from a uniform distribution 
over the square [—10, 10] x [—10, 10]. Table U lists the RMSE of source positions under Gaussian noise, 
with standard deviations cr g aussian = 10~ 3 , 10~ 2 , 10 -1 , and 1. SL£i uses a = 10 6 for the "projector" 
II = (A + <t11 t ) -1 , which is a very conservative value (see Appendix IB1. 



Cgaussian 


SUi 


SLCP 


le-3 


1.5e-3 


1.3e-3 


le-2 


1.41e-2 


le-2 


le-1 


0.1720 


0.1179 


1 


1.7922 


1.4765 



TABLE I: RMSE of SL^i and SLCP under Gaussian noise. 

To compare the algorithms in the presence of outliers, range measurements are created according to 
<k = || • || + Vi, where vi is a Laplacian noise term with standard deviation a^iace £ [0.2, 1.6]. Results 
are also presented for the alternative selective Gaussian outlier generating model of Examples 1-2 with 
^outlier £ [0.5, 2] and <7 gauss i an = 0.04, where outliers only affect measured ranges between the second 
anchor and the source. Tables |IIa] and lllbl list the RMSE of source positions for these two models. We 
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Fig. 3: Mean and uncertainity ellipsoids of MM with different initialization methods and no outliers for 
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Fig. 4: Mean and uncertainity ellipsoids of weighted MM with different initialization methods and outliers, 
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1.8872 
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(a) Laplacian noise (b) Selective Gaussian noise 



TABLE II: RMSE of SLii and SLCP in the presence of outliers. 



conclude from Tables HI and HTl that the relative accuracy of SL£i and SLCP depends on the data generation 
model. For Gaussian noise the RMSE of SL^i is about 30% higher than that of SLCP, whereas in the 
presence of outliers the situation is reversed, and SLCP exhibits an excess RMSE of 10-30%. Interestingly, 
the performance gap is actually larger for selective Gaussian outliers, whose generating model does not 
match the assumptions of SL^i. Similarly to Fig. [2b] it was found that the differences in initialization 
accuracy using SLCP or SL^i are large enough to prevent closing of the RMSE gap after weighted MM 
refinement due to convergence to undesirable extrema of the likelihood function. 

Example 4 [Time Recursive Updating]: This example assesses the performance of the full time- 
recursive procedure, comprising SLCP or SL^i initialization followed by refinement. The network scenario 
has 16 unknown sensors, 4 anchors and 10 target locations, all randomly positioned. A new target 
sighting (the 11th one) becomes available and is processed incrementally, i.e., the position is estimated 
through SLCP or SL^i by fixing all the remaining ones, then all estimates are jointly refined. Results are 
benchmarked against refinement with full batch initialization, which makes a fresh start to the process 
without using any previous knowledge at every new target position to be estimated, solving different and 
increasingly large EDM completion problems for ML initialization. 

This type of incremental approach was used in |[T4l with the SR-LS algorithm of [23] and MM 
refinement for Gaussian noise. SLCP is used here instead of SR-LS because, as shown in 1161 . it increases 
the convergence speed of subsequent iterative methods and also alleviates the problem of convergence 
to local extrema of the ML cost function by providing better initial points than SR-LS does. Figure 
[5] shows the evolution of the Gaussian cost function S1g( x ) during refinement after ranges to the 11th 
target position are sensed (<7 ga ussian = 0.04). The time recursive (SLCP)+MM approach takes advantage 
of previously estimated positions to start with a lower cost than batch (EDM-R)+MM, but reaches the 
same final error value. 
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Behaviour of f2 G (x) with Time Recursive + MM and Batch + MM with a new target position 
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Fig. 5: Evolution of Gaussian cost function ^g( x ) during refinement for EDM-R+MM and SLCP+MM 
approaches, with cr gaLlss i an = 0.04. 

The same network scenario is adopted in the presence of outliers. Figure [6] shows the evolution of cost 
function Ol(x) during refinement for Laplacian outliers (cri ap iacian = 0.1), whose behavior is similar to 
the Gaussian case of Fig. [5] In both Gaussian and Laplacian settings refinement yields similar accuracy 
and convergence speed after batch or time-recursive initializations. Therefore, time-recursive updating is 
seen to retain the essential features of our EDM-based approach to SLAT, namely, very limited need for 
a priori spatial information and fast convergence, at a fraction of the computational cost. 



In this paper, we have presented a ML-based technique to solve a SLAT problem under Gaussian and 
Laplacian noise. An MM method is proposed to iteratively maximize the non-convex likelihood function, 
for which a good initial point is required. Therefore, we have investigated two initialization schemes, 
namely, the batch (EDM) and time -recursive (SLCP/SL^i) approaches that bypass the need for priors on 
sensor/target locations. After the first block of measurements is obtained, an EDM completion method is 
used for the first initialization of the sensor network topology. In our experiments this was accomplished 
reasonably fast (a few seconds) for scenarios with up to about 30 unknown positions. As EDM completion 
is not scalable, we resort to an alternative, lightweight, incremental initialization scheme as additional 
target range measurements become available. The SLCP or SL^i time recursive methods fix the already 



VI. Conclusion 
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Fig. 6: Evolution of Laplacian cost function 0^(x) during refinement for EDM-R-^i+wMM and 
SL^i+wMM approaches, with rjiapiadan = 0.1. 



estimated positions whenever a new position is to be determined; afterwards all positions are given as 
initialization to the likelihood optimization methods. 

Simulation results showed that our method nearly attains the Cramer-Rao lower bound under moderate 
Gaussian noise. In the presence of outliers, EDM-R-^i or SL^i provide more accurate initial position 
estimates than other existing methods. Moreover, when used as input to iterative refinement methods 
they provide a good starting point that reduces the probability of convergence to undesirable extrema, 
yielding improved overall estimation performance. Hence, with this methodology, we obtain a processing 
structure that is robust to outliers and provides a scalable and accurate solution to the SLAT problem. 
Importantly, the algorithms based on t\ norm optimization exhibited this robust behavior in simulation 
not only for Laplacian outliers, but also for an alternative outlier generation technique that did not match 
the underlying Laplacian modeling assumptions. 
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Appendix A 

Convergence of Weighted Majorization Minimization 

To prove (local) convergence of the weighted MM iteration (l22l ) the Laplacian cost function (f2]) is first 
majorized at time t by 

ri(x) = ~ £{4(^(x) - d l3 f + -1} + ~ ^{4iK(x) - 4i) 2 + ^-}, 05) 

i,j V k,j kj 

where fij, g^j and ujj, v\ - are defined in Q and (lTTT ). The inequality ^l(x) < T^(x) follows from 

ri(x) - o L (x) = i ^{4(/^(x) - da) 2 + ^ - 2|^(x) - 



+ 2 J2{ v kj(9kj(x) - d kj ) 2 + ^r~ 2|£?fcj(x) - 4j| } 



(36) 

It is easy to check that ^(x*) = r^(x*), so T^(x) has the properties of a true majorization function for 
the iterate x*. Now the same technique used in (fTOb is applied to majorize (l35l) by a convex quadratic 
function of x, yielding 

^l(x) < \ E{4(/^( x ) - Zdiifijtf) ~ 2dij(V/ -(x*), (x - x*)) + 4) + ^} 
i.j % 3 

1 r H (37) 

+ 9 5j 4(^( X ) - 24iff fci (x*) - 2d J y(Vfly(x t ), (X - X*)) + dy + — . 

As before, equality holds for x = x', so the right-hand side of (l37l) is still a valid majorization function. 
Discarding constant terms the weighted MM iteration (1221 ) results. 

Appendix B 

Properties of Single-Source Localization using SLii 

Proof of Lemma UV-C. 1[ To streamline the notation we define Ki = |||y — bi\\ — and apply the 
KKT condition to the inner optimization problem in d25l ) while fixing y. The Lagrangian function is 

n+l 2 

L(A l7 ) = ^-i + 7 (l T A-l). (38) 

i=i Al 



The KKT conditions are 



^ = -^- + 7* = 0, 1 T A=1. (39) 
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Using d39l ), we find A* = as a solution of the inner optimization problem. Plugging the optimal 

A in the cost function of (T25T ) yields ^Q) 2 , thus establishing the equivalence with (1241 . ■ 
Approximation accuracy of II = lim (T _ >00 (A + all 3 ")" 1 : To decide how large o" should be, let us 
first define n(<r) = (A + all T )- 1 . The norm of the difference to the original definition of II in d29l ) is 
given by 

||n - n(<r)|| F = HA-^K^A- 1 !)- 1 - (1 T A- X 1 + a- x )- x ]l T AT X \\ F 

1 T A- 2 1 



(40) 



(l T A- 1 l)( ( rl r A- 1 l + l)' 
Now assume the most unfavorable case with identical A,- = -Vr, such that 

||n - n(a)|| F = i < e ^ > r^rr - rai- (41) 

cr(n + f) z + 1 (n + tje (n + t) z 

For e = 10 -4 and n + Z = 100, for example, this yields a > 10 2 — 10~ 4 « 10 2 , which is quite low and 
does not raise any numerical issues in commonly available convex optimization solvers. ■ 

Appendix C 
Derivation of CRLB for Gaussian Noise 

The log of the joint conditional pdf for the SLAT problem is (up to an additive constant) 

log/(d|x) = I Yl(\\xi - ej \\ - d i3 f + £(||a fc - ej\\ - d kj f \ , (42) 

[ hi hi ) 

where, similarly to x, d denotes the concatenation of all range measurements. Let us define matrices Mjj 

and Nj that extract individual positions or their differences from the vector of concatenated coordinates 

x as follows 

Mjj-x = Xi — ej, NjX = —ej. (43) 

Thus, (l42l is rewritten as 



log /(d|x) = ~ | X)(||Myx|| - d tj ) 2 + 53(||o fc + N,x|| - d kj f \ 
\ hi hi ) 



(44) 



The first derivative of (|441) with respect to x is 



1 I nt- M^Mj,-x „ Nj(a& + N,x) 

v x io g /(d|x) = - V2 ^(iim^xii - d t ,) " j + + n,x|| - d kj ) r +N \ 

I i,j %3 k,j 3 



(45) 
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The second derivative of (l44l) with respect to x is 




V x log/(d|x) = - < ^ < p^-jp + ^,. x|| I MyMy p-^p 



v- I Nj(q fc + Njx)(afc + N 3 -x) r Nj ||q fe + Njx|| - rf fc j / n t n _ N lK + N jX )(a fc + N,-x) T Nj 
^\ ||a fc + N iX || 2 ' Hafc + NjXll \ j 3 ||a fc + NjX|| 2 



The Fisher information matrix, F x , is obtained by taking the negative expected value of the (l46l ) with 
respect to ranges as lfT3Tl 

tl A\ U 1 [^ M|M^MgMg Nj(a fc +N J x)(a fc + N J xrNj | 

f x = - Ed{ v x log /(d|x) } = -AY, — p-^p — + E kTn^f • 

(47) 

The CRLB matrix in (l34l is taken as the inverse of F x . 
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